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A method is presented for the identification of high-energy neutrinos from gamma ray bursts by means of a large-scale 
neutrino telescope. The procedure makes use of a time profile stacking technique of observed neutrino induced signals 
in correlation with satellite observations. By selecting a rather wide time window, a possible difference between the 
arrival times of the gamma and neutrino signals may also be identified. This might provide insight in the particle 
production processes at the source. By means of a toy model it will be demonstrated that a statistically significant 
signal can be obtained with a km^ scale neutrino telescope on a sample of 500 gamma ray bursts for a signal rate as 
low as 1 detectable neutrino for 3% of the bursts. 
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1 Introduction 

Cosmic radiation is a valuable source of informa- 
tion about various energetic astrophysical processes. 
However, the existence of very energetic cosmic rays 
also raises questions such as : how are they acceler- 
ated and from where do they originate ? 
A variety of possible accelerator mechanisms exists, 
ranging from shock waves produced by exploding 
stars (supernovae) or Gamma Ray Bursts (GRBs) 
to supermassive black holes with strong magnetic 
fields (Active Galactic Nuclei). The current under- 
standing is that protons and electrons are the pri- 
mary particles that are accelerated by electromag- 
netic fields at a cosmic accelerator site. 
In case of a supernova event, a shock is formed 
by the expanding matter envelope when it sweeps 
through the interstellar medium which surrounds 
the exploding star. In such an environment stochas- 
tic processes occur which can accelerate particles to 
very high energies. A detailed treatment [1] shows 
that acceleration by shock waves automatically re- 
sults in a power spectrum, which is in qualitative 



agreement with the observations up to the 'knee' re- 
gion of the cosmic ray spectrum [2] . 
However, this leaves us with the question of which 
events can produce the cosmic rays above the 'knee' 
region. Candidates for the production of the most 
energetic cosmic rays are Active Galactic Nuclei 
(AGN) and Gamma Ray Bursts. The current per- 
ception is that the majority of these objects have a 
similar inner engine, in which infalling matter and 
the likely presence of a strong magnetic field gives 
rise to relativistic shock wave acceleration in two 
back to back jets. 

Interactions of accelerated protons and electrons 
with the ambient photons at the acceleration site 
give rise to very energetic secondary particles, as 
shown in Fig. 1. In particular the pj interactions 
yield a flux of very energetic neutrinos, as depicted 
in more detail in Fig. 2. 

In case of a proton energy of 10^^ eV, i.e. the re- 
gion of the 'knee' of the cosmic ray spectrum, the 
photon energy threshold for A production is about 
10 eV, being the ultraviolet part of the spectrum. 
Since there are many of these UV photons present. 
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Fig. 1. Particle production in jets (courtesy C. Spiering). 
thermal, e synchrotron, inverse Compton 
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Fig. 2. Neutrino production processes. 

the depicted hadronic processes will take place at 
high rates, yielding substantial neutrino fluxes com- 
parable to those of ultrahigh-energy photons. 
In the decay of the A resonance into a nucleon and 
a TT meson, the meson obtains on average 20% of the 
primary proton energy. This yields an average neu- 
trino energy of about 400 TeV for a primary proton 
energy of 10^^ eV. Detailed model calculations [3] 
predict an E~'^ powcrlaw spectrum for the produced 
neutrino flux. Taking into account the fact that the 
atmospheric spectrum is softer [2] and that the neu- 
trino cross section increases with energy [2,4], we 
observe that optimal detection conditions arc ob- 
tained for neutrino telescopes in an energy range of 



about 10-100 TcV [5]. 

Various attempts [6] have been made to identify a 
high-energy neutrino flux in correlation with satel- 
lite observations of GRBs. The performed searches 
for a statistical excess above the background com- 
prise both photon-neutrino coincidence studies and 
investigations of so-called "rolling time windows". 
However, the former will obviously fail in case 
there exists a significant time difference between 
the arrival times of the photon and neutrino fluxes, 
whereas the latter can only be succesful in case 
some GRBs produce multiple neutrino detections 
within the corresponding time windows. So far, no 
positive identifications have been reported. 

From the above it is seen that it would be prefer- 
able to use an analysis procedure that does not re- 
quire the simultaneous arrival of photons and neu- 
trinos and which also provides a high sensitivity in 
case of low signal rates. Such a method, based on a 
time profile stacking technique, is presented here. 



2 The time profile stacking procedure 

In order to obtain a statistical significant result 
even in case of low signal rates, a cumulative proce- 
dure as outlined below has been devised. It is based 
on the generic GRB engine described in the previ- 
ous section, which implies that the arrival times of 
the photons and neutrinos are correlated but are not 
necessarily simultaneous. 

When a GRB is observed by a satellite, the trigger 
time tgrb and burst location on the sky are recorded. 
Afterwards, the data of a neutrino telescope are in- 
spected for a time interval [tgrt — At, tgrb + and 
all arrival times of upgoing muons are recorded rel- 
ative to tgrb- Here At is some predefined time mar- 
gin, which is identical for all observed bursts. An 
upgoing muon is a long reconstructed track in a neu- 
trino telescope pointing backwards to a location in 
the hemisphere opposite to the detector location. 
The usage of upgoing /x tracks allows reduction of 
the (atmospheric) background signals in our analy- 
sis procedure, as outlined lateron. 
For a sample of different GRB observations, the 
above will result in a set of identical time windows 
with upgoing fi arrival time recordings relative to 
the corresponding GRB trigger time. 
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Stacking of all these time profiles will exhibit a 
uniform distribution for background events. How- 
ever, in case the data contain upgoing /i signals cor- 
related with the GRBs, a clustering of data bins is 
expected. Consequently, comparison of the stacked 
time profile contents with a uniform background 
allows the identification of correlated signals. Due 
to the cumulative character of the procedure, large 
statistics can be obtained resulting in a good sensi- 
tivity even in case of low signal rates. 
Any arrival time difference between a photon and 
neutrino signal poses no problem as long as this time 
difference is smaller than At. As such. At should 
be taken as large as possibly allowed by the back- 
ground signals. However, it is obvious that a spread 
in this photon-neutrino arrival time difference will 
reduce the significance of the signal. 

To address the feasibility of the procedure and to 
investigate the effects of the various parameters, a 
toy model which mimics GRB induced signals as 
well as (atmospheric) background has been devised. 
A description of this toy model and the results of 
the above analysis procedure performed on the sim- 
ulated data for a km^-scale detector are presented 
hereafter. 

3 Signal and background generation 

Satellite observations of GRBs have shown [7] 
that the burst locations are homogeneously dis- 
tributed over the sky. Since our analysis procedure 
is based on the detection of upgoing fj, tracks, our 
toy model only generates GRB positions homoge- 
neously distributed over the hemisphere opposite 
to the detector location. 

For each generated burst location we define the 
satellite trigger time to be tgrb = and create a 
time window [—At, At] around it. 
Observations with the AMANDA neutrino tele- 
scope [8] show that a km^-scale detector will ob- 
serve on average 300 upgoing muons per 24 hours 
due to (atmospheric) background, homogeneously 
distributed over the hemisphere. Therefore, each of 
the above time windows will be filled with a number 
of background upgoing muon signals taken from a 
Poissonian distribution with an average number of 

^ See .http://www.phys.uu.nl/~nick/grbmodel 



(300/24) • {2 At /I hour). 

The arrival directions of these background upgoing 
muons are taken to be homogeneously distributed 
over the hemisphere, whereas their arrival times 
are taken to be uniformly distributed in the corre- 
sponding time window. To take the detector time 
resolution at into account, a Gaussian spread with 
a standard deviation at is introduced to the ar- 
rival times. Finally the resulting arrival times are 
recorded as (background) entries in the various 
corresponding time windows. Also the angular po- 
sitions of the arrival directions are recorded, after 
introducing a Gaussian spread corresponding to the 
angular resolution aa of the detector. Recording of 
these angular positions will lateron allow reduction 
of the background by correlation with the actual 
GRB locations. 

By means of a uniform random number generator 
only a fraction / of the generated burst locations is 
selected to yield a single upgoing fj, signal. To mimic 
a time difference r between the photon and neutrino 
burst arrival times, the upgoing /i signal arrival time 
of each signal burst is taken from a Gaussian distri- 
bution with a mean value r and a standard devia- 
tion ar- Before these signal muon arrival times are 
added to the corresponding time windows, a Gaus- 
sian spread corresponding to the detector time res- 
olution at is introduced. The arrival directions of 
these signal upgoing muons are recorded as the lo- 
cations of the corresponding bursts, after introduc- 
ing a Gaussian spread corresponding to the detector 
angular resolution aa- 

Introduction of realistic values for the various toy 
model parameters outlined above will allow to inves- 
tigate the feasibility of detecting neutrino induced 
signals in a large scale neutrino telescope. In an ac- 
tual experimental data analysis effort one obviously 
has to account for several additional (systematic) ef- 
fects like detector stability, track reconstruction ef- 
ficiency and so on. However, these are detector spe- 
cific effects and fall beyond the scope of the present 
studies. 



4 Analysis of simulated data 

The only large scale neutrino telescope currently 
in operation is IceCube [9] and as such we use the 
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parameters of this detector [5] as benchmark values 

for our present studies. 

The expected data rates for the full km'^ scale detec- 
tor allow a time margin At of 1 hour. This implies 
an average number of background signals of about 
25 upgoing muons for each individual GRB time 
window. The time resolution for the reconstructed 
muon tracks will be of the order of the time it takes 
for a muon to cross the detector volume. As such 
wc take at = 10 /is as a conservative estimate for 
the detector time resolution. Experience with the 
analysis of the AMANDA data [8,11] together with 
detector simulation studies [5] show that a realistic 
estimate for the angular resolution is obtained by 
taking Ua = 1°. 

The remaining parameters of our toy model are re- 
lated to the characteristics of the various bursts. 
Based on the processes sketched in Fig. 1, a reason- 
able estimate for the possible photon-neutrino ar- 
rival time difference and its spread can be obtained 
from the actual burst duration. Satellite observa- 
tions [7] exhibit a mean burst duration of about 
30 seconds. As such we take r = 30 s and ct,- = 30 s. 
As mentioned before, for evaluation of the currently 
presented procedure the value of r is actually irrel- 
evant as long as it is smaller than the time margin 
At. 

This leaves us with only two free parameters : the 
fraction / of GRBs that actually induces an upgoing 
muon signal and the bin size to be used for the time 
profiles. In order to optimise the time bin clustering 
of the signals, the bin size should be taken to be of 
the order of the temporal signal spread cxr . However, 
since the observed redshifts of GRBs [7] exhibit a 
median value of z = 1.9 with a spread of 1.3, cos- 
mological time dilation effects have to be taken into 
account. It should be noted, however, that in case 
both the photon and neutrino production processes 
are taking place continuously throughout the jet ex- 
istence, the cosmological time dilation is already in- 
cluded in the observed gamma burst duration and 
consequently also in the temporal signal spread a^. 
Nevertheless, we always account for a possible ad- 
ditional cosmological time dilation and take for the 
time profile bin size a conservative value of 5(7,-, cor- 
responding to 150 s. 

It should be noted here that restricting the analysis 
to short duration bursts allows for smaller time bins 



and consequently more detailed time profile studies. 
The fraction / we keep as a free parameter in order 
to determine the sensitivity of our analysis proce- 
dure for different sizes of the GRB sample. 

For a first investigation of the performance of 
the procedure we generated 100 GRBs in one hemi- 
sphere. This corresponds to about 2 years of opera- 
tion of the Swift satelHte [10], which currently is the 
main source of GRB triggers. All parameters were 
set to the values mentioned above and for the frac- 
tion / we used a value of 10% [12]. The resulting 
stacked time profile is shown in Fig. 3. 
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Fig. 3. Stacked time profile for 100 GRBs with / = 0.1. 
Further details can be found in the text. 

Since in our toy model we have access to all in- 
formation, we are also able to construct the corre- 
sponding stacked time profile from the background 
signals only. This background stacked time profile 
is shown in Fig. 4. 
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Fig. 4. Stacked time profile corresponding to the back- 
ground data of Fig. 3. 
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Comparison of the number of entries from Fig. 3 
and Fig. 4 shows that 8 of our generated GRBs in- 
duced a signal in the stacked time window. How- 
ever, due to the presence of a large background we 
are not able to identify the GRB signals on the basis 
of our observations of Fig. 3 alone. 

Reduction of the background without significant 
signal loss can be achieved by only investigating a 
certain angular region centered around the actual 
GRB position. As detector angular resolution we 
have (Ja = 1°, so restricting ourselves to an angular 
region of 5° around the GRB location will reduce 
significantly the background while preserving basi- 
cally all signal muons. 

The stacked time profile of our previous generation, 
but now restricted to an angular region of 5° around 
the burst location, is shown in Fig. 5. 




-3000 -2000 -1000 1000 2000 3000 

Upgoing jj arrival time (in sec) w.r.t. GRB y trigger 

Fig. 5. Stacked time profile for 100 GRBs with / = 0.1 
and restricted to an angular region of 5° around the 
actual burst location. 



Visual inspection of Fig. 5 raises some doubts to 
a conclusion that the observed time profile results 
from a uniform background distribution. This is 
confirmed if we investigate the corresponding back- 
ground distribution as shown in Fig. 6. 

Comparison of Fig. 5 and Fig. 6 allows the iden- 
tification of the GRB signals in the central bin. In 
the analysis of experimental data, however, we don't 
have access to the actual corresponding background 
distribution. As such, we need to quantify our de- 
gree of (dis)belief in a background observation solely 
based on the actually recorded signals like in Fig. 5. 
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Fig. 6. Stacked time profile corresponding to the back- 
ground data of Fig. 5. 

5 Bayesian assessment of the significance 

Consider two propositions A and B and some 
prior information /. We introduce the notation 
p{A\BI) to represent the probability that A is true 
under the condition that both B and / are true. 
Following the arguments of extended logic [13] we 
automatically arrive at the so-called theorem of 
Bayes 



p{B\AI) = p{B\I) 



p{A\BI) 

Pirn 



(1) 



The above theorem is extremely powerful in the pro- 
cess of hypothesis testing, as will be shown here. 
Consider a hypothesis H in the light of some ob- 
served data D and prior information /. By we 
denote an unspecified alternative to H. This implies 
that H^: is just the proposition that H is false. From 
eq. (1) we immediately obtain 



p{H\DI) _ p{H\I) p{D\HI) 
p{H,\DI) ~ p{H,\I) p{D\HJ) 



(2) 



Introducing an intuitive decibel scale, we can ex- 
press the evidence e{H\DI) for H relative to any 
alternative based on the data D and prior informa- 
tion / as : 



e{H\DI) = lOlog^o 



p{H\DI) 



[p{H,\DI)_ 
Combined with eq. (2) this yields 

e{H\DI) = e{H\I) + lOlogio ^^^^'^^^ 



p{D\HJ) 



(3) 



(4) 



5 



To quantify the degree to which the data sup- 
port a certain hypothesis H, we introduce the 
Bayesian observables il) = — 101og;^oP(-^l-^-^) 
V"* = — 10 logiQp{D\H^:I). Since the value of a prob- 
abihty always lies between and 1, we have ip ^ 
and = 0. Together with eq. (4) we obtain 

e{H^\DI) = e{H^\I)+il^-'ilj^^e{H^\I)+^ . (5) 

In other words : there is no alternative to a certain 
hypothesis H which can be supported by the data 
D by more than ip decibel, relative to H. 
So, the value i/j = — lOlogiQ p{D\ HI) provides the 
reference to quantify our degree of belief in H. 

In our evaluation of the stacked time profile the 
main question is to which degree we believe our 
observed distribution to be inconsistent with re- 
spect to a uniform background. This question can 
be answered unambiguously if we are able to deter- 
mine the ip value corresponding to the uniform back- 
ground hypothesis based on our observed stacked 
time profile. 

The process of recording background signals is 
identical to performing an experiment with m dif- 
ferent possible outcomes {^i, Am} at each trial. 
Obviously, m is in our case just the number of bins 
in the time profile and the number of trials n is the 
number of entries. 

In case all the probabilities pk corresponding to the 
various outcomes on successive trials are inde- 
pendent and stationery, the experiment is said to 
belong to the Bernoulli class Bm [13]. It is clear that 
our data recordings according to a uniform back- 
ground hypothesis satisfy the requirements of B^- 
The probability p(ni . . .n„i\B„J) of observing 
occurrences of each outcome A^ after n trials is 
therefore given by the multinomial distribution [13] . 
Consequently, the probability for observing a spe- 
cific set of background data D consisting of n entries 
is given by 



p{D\BmI) 



I Pi Pm 



(6) 



ni! ■■■rin. 

This immediately yields the following expression for 
the i/j value according to a uniform background hy- 
pothesis 



ijj = -10 



logion! + ^(nfclogioPfe - logio ^k^-) 

k=l 



When a signal from a uniform background is being 
recorded in our time window, there is no preference 
for any specific time bin. This implies that in our 
case all pk values are identical and equal to ■m~^. As 
such we can evaluate the '0 value of eq. (7) for any 
set of observed data D. 

5. 1 Relation to a frequentist approach 

In the case of large statistics we can use Stirling's 
approximation In x! = xlnx — x for x 3> 1 in eq. (7) . 
Together with the fact that J2f^k = n this yields the 
frequentist approximation 



wJ2^klog^J — ) ■ (8) 



Furthermore, for a "near match" scenario we have 
Hk ~ npk- In such a case we can use the series ex- 
pansion In a; = (x — 1) — ^ (x — 1)^ -|- . . . , which yields 



^nfeln 



k=l 



f nk \ 
\npkJ 



-I m / 

1 [nk 



npk) 



k=l 



npk 



(9) 



This yields the correspondence with the statistic 

{nk - npkY 



k=i 



npk 



(10) 



(7) 



Equation (10) allows a frequentist evaluation 
of the statistical significance of our observations. 
However, this will only provide meaningful results in 
case the conditions mentioned above are satisfied. In 
case a rather unlikely event happens to be observed 
within a small number of trials, a analysis may 
lead to completely wrong conclusions whereas the 
Bayesian approach outlined above will provide the 
correct results [13]. As such, the present studies will 
be based on the exact Bayesian expression of eq. (7) . 

6 Discovery potential 

Evaluation of the expression of eq. (7) for the data 
displayed in Fig. 3 yields ip = 713.38 dB. Since these 
data don't allow the identification of a GRB signal, 
this rather high ip value must be due to background 
fluctuations. This is indeed confirmed by investiga- 
tion of the corresponding background data shown in 
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Fig. 4, which yield iptkg = 709.43 dB. Consequently, 
it is required to determine the ip value of the cor- 
responding background before the statistical signif- 
icance of an observed time profile can be evaluated. 

In our toy model studies we have directly access 
to the corresponding background time profile, but 
this will in general not be the case in an actual 
experimental data analysis effort. One way to inves- 
tigate background signals is to record data as out- 
lined above, but with fictativc GRB trigger times 
not coinciding with the actual tgrb- This method we 
call "on source off time". In order to have similar 
detector conditions for both the signal and back- 
ground studies, the fictative trigger times should 
be chosen not too distinct from the actual tgrb- 
Recording background data in a time span covering 
1 day before and 1 day after the GRB observation 
will allow the investigation of at least 25 different 
background time profiles per burst. These in turn 
will yield the corresponding different stacked back- 
ground time profiles which allow the determination 
of an average value i/Jbkg and the corrresponding 
root mean square deviation Sbkg- 

The processing of extra background data as de- 
scribed above might turn out to become unprac- 
tical, due to e.g. data volume. In such a case one 
might consider using the remaining data of the off 
source locations of the actual time windows. Such a 
method we call "off source on time" . 
The performance of the "off source on time" method 
obviously depends on various detector conditions, 
which may limit the feasibility of such a background 
determination. To overcome these possible limita- 
tions, one could also envisage using the actual ob- 
served time profile and randomise the entries in 
time. By performing several randomisations, a rep- 
resentation of the corresponding background is ob- 
tained. This method we call "time shuffling". It 
should be noted, however, that in the case of a large 
signal contribution the time shuffling method will 
underestimate the significance of the signal. 

In view of the above, we will use the "on source 
off time" method in our toy model studies by gen- 
erating 25 different background samples and per- 
forming our analysis procedure for each of them. 
In the case of the situation reflected by Fig. 3 this 
yields 'ipbkg = 692.04 dB and Sbkg = 21.19 dB, which 
is seen to be in excellent agreement with the actual 



background value corresponding to Fig. 4. 
Comparison of the actually observed value of 
713.38 dB with the reconstructed background val- 
ues immediately shows that no significant signal is 
observed. 

However, evaluation of the data corresponding to 
Fig. 5 yields tp = 218.78 dB with background val- 
ues -ijjbkg = 99.62 dB and Sbkg = 23.98 dB. Here a 
statistically significant signal is obtained. 

For a uniform background and large statistics, the 
Bayesian tp observable can be approximated by the 
frequentist statistic, as indicated in eqs. (8)-(10). 
This implies that the statistical significance for de- 
viation from a uniform background distribution can 
be expressed in terms of a standard deviation a by 
comparison of the actually observed i/j value of the 
stacked time profile with the corresponding 'ipbkg and 
Sbkg background values. This is illustrated in Fig. 7 
for a sample of 250 different background samples 
according to the situation reflected in Fig. 5. The 
distribution of the obtained ipbkg values as shown in 
Fig. 7 exhibits a Gaussian profile with a mean value 
and standard deviation which are indeed consistent 
with the above ipbkg and Sbkg values, respectively. 
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Fig. 7. Distribution of ipbkg values for a large sample 
of background distributions according to the situation 
depicted in Fig. 5. 

Variation of the number of GRBs allows a deter- 
mination of the minimal value of the fraction / for 
which a statistically significant signal can be ob- 
tained. Common practice is to claim a discovery in 
the case a significance in excess of 5a is obtained. 
Following the procedure outlined above this leads 
to the discovery sensitivities as shown in Fig. 8. 
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Fig. 8. Sensitivities corresponding to a 5a signal signifi- 
cance. 

It should be noted that the actually achievable 
sensitivities are depending on various detector spe- 
cific parameters and the quality of the available 
data. As such, all parameters as well as the amount 
of possible background samples will have to be op- 
timised for each specific experimental data analysis 
scenario. 

In case no significant signal can be identified from 
an experimentally observed stacked time profile, 
values like the ones presented in Fig. 8 provide the 
basis for a fluence limit determination. 



7 Summary 

The method introduced in this report allows iden- 
tification of high-energy neutrinos from gamma ray 
bursts with large scale neutrino telescopes. The pro- 
cedure is based on a time profile stacking technique, 
which provides statistical significant results even in 
the case of low signal rates. 

The performance of the method has been inves- 
tigated by means of toy model studies based on 
realistic parameters for the future IceCube km^ 
neutrino telescope and a variety of burst samples. 
From these investigations it is seen that a 5cr signif- 
icance is obtained on a sample of 500 bursts with a 
signal rate as low as 1 detectable neutrino for 3% 
of the bursts. 

Finally, it should be realised that the actually 
achievable sensitivities are depending on various 
detector specific parameters and the quality of the 
available data. These aspects, however, fall beyond 



the scope of the present report. 
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